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Abstract. Consideration is given to the influence of an underwater landslide on 
waves at the surface of a shallow body of fluid. The equations of motion which 
govern the evolution of the barycenter of the landslide mass include various dis- 
sipative effects due to bottom friction, internal energy dissipation, and viscous 
drag. The surface waves are studied in the Boussinesq scaling, with time-dependent 
bathymetry. A numerical model for the Boussinesq equations is introduced which 
is able to handle time-dependent bottom topography, and the equations of motion 
for the landslide and surface waves are solved simultaneously. 

The numerical solver for the Boussinesq equations can also be restricted to imple- 
ment a shallow-water solver, and the shallow-water and Boussinesq configurations 
are compared. A particular bathymetry is chosen to illustrate the general method, 
and it is found that the Boussinesq system predicts larger wave run-up than the 
shallow- water theory in the example treated in this paper. It is also found that the 
finite fluid domain has a significant impact on the behavior of the wave run-up. 
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1. Introduction 

Surface waves originating from sudden perturbations of the bottom topography 
are often termed tsunamis. Two distinct generation mechanisms of a tsunami are 
underwater earthquakes, and submarine mass failures. Among the broad class of 
submarine mass failures, landslides can be characterised as translational failures which 
travel considerable distances along the bottom profile [30, 45] . In the past, the role of 
landslides and rock falls in the excitation of tsunamis may have been underestimated, 
as most known occurrences of tsunamis were accredited to seismic activity. However, 
it is now more accepted that submarine mass failures also contribute to a large portion 
of tsunamis [51], and recent years have seen a multitude of works devoted to the study 
of such underwater landslides and the resulting effect on surface waves [3, 12, 14, 21, 
29, 30, 39, 40, 44, 51]. As suggested in [23], it is possible for underwater landslides 
and earthquakes to act in tandem, and produce very large surface waves 
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A natural question to ask is whether the effect of underwater landshdes on surface 
waves can be such that they may pose a danger for civil engineering structures located 
near the shore. Consequently, one important issue is the wave action and in particular 
the run-up and draw-down at beaches in the vicinity of the landslide. While the draw- 
down itself may not pose a threat, one consequence of a large draw-down can be the 
amplification of the run-up of the following positive wave crest [17, 49]. 

There have been many numerical and a few experimental studies devoted to this 
subject, but it is generally difficult to include many of the complex parameters and 
dependencies of a realistic landslide into a physical model. Therefore, most workers 
attempt to distill the problem to a model setup where many effects such as turbulence 
and sedimentation are disregarded. For example, Grilli and Watts [30] study tsunami 
sensitivity to several landslide parameters in the case of a landslide in a coastal area of 
an open ocean. In particular, dependence on the landslide shape and the initial depth 
of the landslide location are studied, and it is found that the landslide with the small- 
est length produced the largest waveheight and run-up, and that the wave run-up at 
an adjacent beach is inversely proportional to the initial depth. The work in [30] relies 
on integrating the full water-wave equations using an irrotational boundary-element 
code, and using an open boundary with transmission conditions [28, 27]. While most 
works have considered a given dynamics for the landslide, the bottom motion in [30] 
is described by an ordinary differential equation similar to the one used here. Thus 
the motion of the landslide is computed using a differential equation derived from 
first principles using Newtonian mechanics. However to expedite comparison with 
experiments, the landslide in [30] is considered moving on a straight inclined bottom 
with constant slope. 

More recently, Khakimzyanov and Shokina [35], and Chubarov et. al. [12] have 
also used a differential equation to find the bottom motion. One major novelty in 
their work is that the landslide motion is computed on a bottom with an arbitrary 
shape. The time-dependent bathymetry is then used to drive a numerical solver of the 
shallow-water equations. An advantage of this approach when compared to [30] is the 
reduced computation time. On the other hand, the description of the wave motion in 
the shallow-water theory is only approximate, and in particular, one important effect 
of surface waves, namely the influence of frequency dispersion is neglected. 

The main aim of the current work is to introduce the effect of dispersion into the 
description of the wave motion at the surface of the fluid while keeping the simplicity 
of the shallow-water approach. To this end, we utilise the so-called Peregrine system 
which is a particular case of a general class of model systems which arise in the 
Boussinesq scaling [10]. A common feature of all Boussinesq-type systems is that 
they allow a simplifled study of surface waves in which both nonlinear and dispersive 
effects are taken into account. In the present case, we need to use a Boussinesq system 
which can handle complex and time-dependent bottom topography. Such a system 
was derived by Wu [56] , and can be used in connection with the dynamic bathymetry. 

We conduct two main experiments. First, a comparison with the shallow- water 
theory is carried out. Second, the dependence of the tsunami characteristics on the 
initial depth of the landslide is investigated. The main flndings of the present work 
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Figure 1. This figure shows the fixed bathymetry z - ho^x) and the po- 
sition of the landslide after 50 s. The position of the barycenter is indicated 
by a black dot. 

are that the predictions of the shallow-water and Boussinesq theory are divergent for 
the cases treated in this paper, and that the effect of a finite fluid domain, such as a 
river, lake or fjord [44] can lead to significantly different behavior when compared to 
tsunamis on an open ocean. 

The Boussinesq model in this paper is based on the assumption of an inviscid 
fluid, and irrotational flow. These are standard assumptions in the study of surface 
waves, and generally give good results, unless there are strong background currents 
in the fluid. Another effect which is not taken account of here is the wave resistance 
on the landslide due to waves created by the motion of the landslide. However, 
as observed in [32], this effect is negligible for most realistic cases of underwater 
landslides. Viscosity is included in the dynamic model for the landslide as will be 
shown in the next section. In order to capture the effect of slide deformation during 
the evolution, a damping term in the equation of motion is included to model the 
internal friction in the landslide mass. 

The paper is organised in the following way. In Section 2, the equation of motion 
for the landslide is developed. Then in Section 3, the Boussinesq model is recalled. 
In Section 4, solitary-wave solutions of the Peregrine system are found numerically. 
In Section 5, the numerical scheme for the Boussinesq system is explained and the 
numerical method is tested using the exact solutions of Section 4. Section 6 contains 
results of numerical runs for a few speciflc cases of bottom bathymetry, a parameter 
study of wave run-up in relation to the initial depth of the landslide, and a comparison 
with the shallow-water theory. 



In this section we briefly present a mathematical model of underwater landslide 
motion. This process has to be addressed carefully since it determines the subsequent 
formation of water waves at the free surface. In the present study, we will assume 



2. The landslide model 
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the movable mass to be a solid body with a prescribed shape and known physical 
properties. Since the landslide mass and volume is preserved during the evolution, 
it is sufficient to determine the position of the barycenter x - xdt) as a proxy for 
the motion of the whole body. As observed in the introduction, most studies of wave 
generation due to underwater landslides are based on prescribed bottom motion, or on 
solving the equation of motion on a uniform slope while taking account of different 
types of friction and viscous terms. Examples of such works are [13, 42, 54]. A more 
general approach was recently pioneered by Khakimzyanov and Shokina in [35], where 
curvature effects of the bottom topography were taken into account. Since this model 
is applicable to a wider range of cases, we follow the approach of [35]. 

The static bathymetry is prescribed by a sufficiently smooth single- valued function 
z = -hQ(x), and the landslide shape is initially prescribed by a localised function 
z = Co (a;). To be specific, in this study we choose the following shape function for the 
landslide mass: 

Ux) = a\ k(^^-osC-^)), \x-xo\<l ^2.1) 
[ 0, \x- xo\ > |. 

In this formula, A is the maximum height, i is the length of the slide and xq is the 
initial position of its barycenter. It is clear that the model description given below 
and the method of numerical integration used in the present work is applicable to 
any other reasonable shape. 

Since the landslide motion is translational, its shape at time t is given by the 
function z = C{x,t) = (oi^ ~ Xc{t)). Recall that the landslide center is located at a 
point with abscissa x = Xc{t). Then, the impermeable bottom for the water wave 
problem can be easily determined at any time by simply superposing the static and 
dynamic components. Thus the bottom boundary conditions for the fluid are to be 
imposed at 

z = -h(x, t) = -ho{x) + ({x, t). 

To simplify the subsequent presentation, we introduce the classical arc-length pa- 
rameterisation, where the parameter s - s{x) is given by the formula 

s = L{x)= r^l + iKiOydt (2.2) 

Jxo 

The function L{x) is monotone and can be efficiently inverted to yield the original 
Cartesian abscissa x = L'^(s). Within the parameterisation (2.2), the center of the 
landslide is initially located at a point with the curvilinear coordinate s = 0. The 
local tangential direction is denoted by r and the normal direction by n. 

A straightforward application of Newton's second law reveals that the landslide 
motion is governed by the differential equation 

d^s 

where m is the landslide mass and i^rCO is the tangential component of the sum of 
forces acting on the moving submerged body. In order to project the forces onto the 
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axes of the local coordinate system, the angle 9{x) between r and Ox is needed. This 
angle is determined by 

6(x) = arctan(/iQ(x)). 

Let us denote by and pi the densities of the water and landslide material cor- 
respondingly. If V is the volume of the slide, then the total mass m is given by the 
expression 

m:= {pi + c^,p^)V, (2.3) 

where is the added mass coefficient. As explained in [6], a portion of the water mass 
has to be added to the mass of the landslide since it is entrained by the underwater 
body motion. For a cylinder, the coefficient is equal exactly to one, but in the 
present case, the coefficient has to be estimated. The volume of the sliding material 
is given hj V = W ■ S, where W is the landslide width in the transverse direction, and 
5* can be computed by 

S= (o(x)dx. 

JR 

The last integral can be computed exactly for the particular choice (2.1) of the land- 
slide shape to give 

V = -£AW. 
2 

The total projected force acting on the landslide can be conventionally represented 
as a sum of the force Fg representing the joint action of gravity and buoyancy, and 
the total contribution of various dissipative forces. 

The gravity and buoyancy forces act in opposite directions and their horizontal 
projection Fg can be easily computed by 

Fg{t) = (Pi - Pw)Wg C{x, t) sin(^(x)) dx. 

Now, let us specify the dissipative forces. The water resistance to the motion of 
the landslide Fj. due to viscous dissipation is proportional to the maximal transverse 
section of the moving body and to the square of its velocity. In addition, the coeffi- 
cient sign(^) is needed to dissipate the landslide kinetic energy independently of its 
direction of motion. Thus the force Fr takes the form 



where Cd is the resistance coefficient of the water. The friction force Ff is proportional 
to the normal force exerted on the body due to the weight: 

Ff^-Cf sign|^jAr(x,t). 

The normal force N{x, t) is composed of the normal components of gravity and buoy- 
ancy forces, but also of the centripetal force due to the variation of the bottom slope: 

N{x,t) = {pi - Pw)gW [ C{x,t)cos[9{x)) dx + peW f C{x,t)K{x)i'^\ dx. 
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Figure 2. This figure shows position and velocity of the barycenter of the 
landslide as functions of dimensional time for three different values of the 
friction coefficient c/. 



Here k{x) is the signed curvature of the bottom which can be computed using the 
formula 



k(x) 



(i + {Kix)y) 



3 ■ 

2 



We note that the last term vanishes for a plane bottom since k{x) e in this particular 
case. Energy loss inside the sliding material due to internal friction is modeled by 

ds 

F,^-c,p,WS-, 
at 

where is an internal friction coefficient. Finally, dissipation in the boundary layer 
between the landslide and the solid bottom is taken account of by the term 



ds 



dt 



where is the Chezy coefficient. Finally, if we sum up the contributions of all the 
forces described above, we obtain the second order differential equation 



(7 + c^)S^ = (7 - i)9(h(t) - Cfam(t)) 



dt^ 



/ x/' ^ / X 1 .\/'ds\2 ^ds „ds 
ait){cnh{t) ^ -c,A){-) -c^^S--c,i- 



ds 

di 



(2.4) 
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Figure 3. This figure shows the velocity and acceleration of the barycenter 
of the landslide as a function of dimensional time. The friction coefficient is 
Cf = tan(3°). The discontinuities in the acceleration are due to the coefficient 
sign(^) in the definition of the friction force. 

where 7 := ^ > 1 is the ratio of densities, a{t) := sign^^^ and the integrals /i, 2,3(0 
are defined by 

Ii{t)= f C{x,t)sm(9{x))dx, 
h(t) = / C(x,t)K(x) dx. 

JR 

In order to obtain a well-posed initial value problem, equation (2.4) has to be supple- 
mented with initial conditions for s(0) and s'(0). In the remainder we always take 
homogeneous initial conditions, and consider the motion driven only by the gravi- 
tational acceleration of the landslide. However, different boundary conditions might 
also be reasonable from a modeling point of view. 

In order to approximate solutions of equation (2.4), we employ the Bogacki-Shampine 
third-order Runge-Kutta scheme. The integrals /i,2,3(i) are computed using the trape- 
zoidal rule, and once the landslide trajectory s = s{t) is found, we use equation (2.2) 
to find its motion x = x(t) in the initial Cartesian coordinate system. This yields the 
bottom motion that drives the fluid solver. 

3. The Boussinesq model 

Once the motion of the landslide is determined, and therefore the time-dependent 
bathymetry h{x,t) = ho{x) -C{x,t) is given, the next step is to consider the coupling 
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between the bathymetry variations and the evolution of surface waves. The main 
assumptions on the fluid are that it is inviscid and incompressible, and that the 
flow is irrotational. Under these assumptions, the potential-flow free surface problem 
governs the motion of the fluid. However, in the present case, the fluid is shallow, 
and the waves at the surface are of small amplitude when compared to the depth of 
the fluid. In that case, the potential-flow problem may be simplified, and the model 
used in this paper is a variant of the so-called classical Boussinesq system derived by 
Boussinesq [10]. 

Let us first consider the case of an even bottom, and a constant fluid depth do. 
Denote a typical wave amplitude by a, and a typical wavelength by A. The parameter 
a = ^ then describes the relative amplitude of the waves, and the parameter = p- 
measures the 'shallowness' of the fluid in comparison to the wavelength. In the case 
when both a and /3 are small and approximately of the same order of magnitude, the 
system 

r]t + doUx + {rju)^ = 0, 

d^ (3.1) 

ut + gr]x + uu.^ - Y'lJ'xxt = 

may be used as an approximate model for the description of the evolution of the 
surface waves and the fluid flow. In (3.1), t] denotes the deflection of the free surface 
from its rest position, and u denotes the horizontal fluid velocity at a height z = 
do(-l + 1^1/3) in the fluid column if z is measured from the rest position of the free 
surface. The same equation appears if the velocity is taken to be the average of the 
horizontal velocity over the flow depth. 

The system (3.1) was first derived by Peregrine in [43], and falls into a general 
class of Boussinesq systems, as shown in the systematic studies [9, 38]. As opposed 
to the shallow-water approximation, the pressure is not assumed to be hydrostatic, 
and the horizontal velocity varies with depth. In fact, the horizontal velocity profile 
is a quadratic function of z [55]. 

The derivation of (3.1) given in [43] also featured and extension to non-constant 
but time- independent bathymetry. However, the present case of a dynamic bottom 
profile calls for a system which allows for time-dependent bathymetry, and such a 
system was derived in [56]. Given a bottom topography described hj z - -h{x,t), 
the system takes the form 

Vt + {hu)x + (?7m)x + /it = 0, 

Ut + gr}x + uu^ = —u^xt + -h {ht + {hu)^)^^ . 

In order for this system to be asymptotically valid, we need a ~ /3 as before. Moreover, 
concerning the unsteady bottom profile, we make the assumptions that hx < 0(a/3^/^), 
and ht<0(a/3y^). 

In comparison to the shallow-water equations with a time-dependent bottom to- 
pography, the system (3.2) has additional terms on the right-hand side of the second 
equation. The effect of these terms is to incorporate frequency dispersion into the 
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model. One practical aspect of this modification is that wave breaking can be com- 
pletely avoided as long as the amplitude of the waves is small enough. Wave breaking 
is also possible in evolution systems of Boussinesq type [7], but the amplitudes oc- 
curring in the present problem are far from the breaking limit. One may argue that 
the dispersion in (3.2) is too strong in comparison with dispersion in realistic water 
waves. However, the linear dispersion relation of (3.2) with fixed even bottom is 
closer to the dispersion relation of the original water-wave problem than most other 
standard Boussinesq equations [7]. 

4. Solitary waves 

Before the numerical method for approximating solutions of (3.2) is presented, we 
digress for a moment, and explain how to find numerically exact solutions of the 
system (3.1). These solutions will later be used to test the implementation of the 
numerical procedure. Assuming the special form 

and substituting this representation into the governing equations (3.1), there appears 

-Csf]' + ({d + r])u)' = 0, 

1 

-Csu' + -{u^ + gr]' + CsY^"' = 0. 

Assuming decay of both r] and u to zero as oo, the integration of the mass 

conservation equation from -oo to ^ gives the following relation between t] and u: 

u = -^, V = . (4.1) 

a + 7] Cg- u 

The momentum balance equation can now be integrated to yield 

/ 2 

Finally, in order to obtain a closed form equation in terms of the velocity u, we 
substitute the expression (4.1) for r] into (4.2). The resulting differential equation 
can be written in operator notation as 

Cu=Af{u), 

where the linear operator C and the nonlinear operator Af, are defined respectively 
by 

//\ J AT/ ^ 1 2 . gdu 



,(^u--u") + -u' + grj = 0. (4.2) 



Cu := Cs(u - —u"\ , and Af(u) ■■= -u^ + 



3 ' 2 Cg - u 

While nothing formal appears to be known about existence of localised solutions of 
(4.1) (4.2), it is straightforward to compute approximations of solitary waves numer- 
ically. In particular, one may use the well known Petviashvili iteration method which 
takes the form 

.M{u^)-\ ^^^^^^^^ I . (4.3) 



10 



D. DUTYKH AND H. KALISCH 




Figure 4. This figure shows the comparison of the numerical approxima- 
tion of solitary wave solutions of (3.1) to Grimshaw's third-order asymptotic 
approximation of solitary waves using the Euler equations for the full water 
wave problem. The left panel shows the surface elevation, and the right 
panel shows the horizontal velocity at z - do{~l + a/i/3). 



The exponent q is usually defined as a function of the degree p of the nonlinearity, 
with the rule of thumb that the expression q := ^ generally works well. In our case, 
the nonlinearities are quadratic, so that we choose p = 2, and hence q - 1. 

The Petviashvili method was analyzed in [41], and can be very efficiently imple- 
mented using the Fast Fourier Transform [22]. The iteration can be started for in- 
stance with the third-order asymptotic solution of Grimshaw [31]. The iterative 
procedure is continued until the Loo norm between two successive iteration is on the 
order of machine precision. Figure 4 shows approximate solitary-wave solutions of 
(3.1) with various wave speeds, and compares them to the third-order asymptotic 
approximation of solitary-wave solutions of the full water-wave problem obtained by 
Grimshaw [31]. The left panel shows comparisons of the free-surface excursion, while 
the right panel shows a comparison of the horizontal component of the velocity field, 
evaluated at the non-dimensional height z given by 5 = -1 + \/Tj3. Figure 5 shows a 
comparison of the wavespeed-amplitude relation between the solitary-wave approxi- 
mation of (3.1) and the ninth-order asymptotic approximation to the full water-wave 
problem obtained by Fenton [20]. 



5. The numerical scheme 

For the numerical discretisation, a finite-volume discretisation procedure similar to 
the one used in [4, 5] is employed. Let us take as a unit of length the undisturbed 
depth do of the fluid above the barycenter of the landslide, and as a unit of time the 
ratio \J^- Then the Peregrine system (3.2) is rewritten in terms of the total water 
depth H as 

Ht + [Hu], = 0, (5.1) 

Ut + [W + {H - h)] = \hh^u + \h{hu)^^t - l:h'^u^xt, (5.2) 
•^2 2 b 
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Figure 5. This figure shows the amphtude-speed relation of soHtary wave 
solutions of (3.1) and of Fenton's ninth-order asymptotic approximation of 
solitary waves using the Euler equations for the full water wave problem. 

The system (5.1), (5. 2) can be formally rewritten in the form 

+ [F(V)], = §b + M(V), 
where the density V and the advective flux F(V) are deflned by 



(5.3) 



V 



F(V) 



Hu 



The source term is defined by 



^■^hhxtt I ' 



and the dispersive term is defined by 

M(V) : 







We begin our presentation by a discretisation of the hyperbolic part of (5.1), (5.2), 
which is the classical nonlinear shallow-water system, and then discuss the treatment 
of dispersive terms. The Jacobian of the advective fiux F(V) is easily computed to 
be 

d¥{Y) _ /u H\ 



A(V) 



dV \^ 

and it is clear that A(V) has the two distinct eigenvalues 



M ± Co 



H. 



The corresponding right and left eigenvectors are the columns of the matrices 



R 



'H 



L = R-^ 



1 
2 



-1 



H 



c-\ 
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We consider a partition of the real line M into cells (or finite volumes) Ci = [a^j-i , 
with cell centers Xi = + x-^i) {i eZ). Let Axi denote the length of the cell Ci. 

In the sequel we will consider only uniform partitions with Axi = Ax, Vi 6 Z. We 
would like to approximate the solution Y{x,t) by discrete values. In order to do so, 
we introduce the cell average of V on the cell Ci (denoted with an overbar), i.e., 

Viit) = {H,{t),u,{t)) = ^ f Y{x,t)dx. 

A simple integration of (5.3) over the cell Ci leads to the exact relation 

^ + ^[F(V(x,,.,t))-F(V(x,_.,t))] = l-£s,(Y)dx . 

Since the discrete solution is discontinuous at cell interfaces Xj^i (z e Z), we replace 
the flux at the cell faces by the so-called numerical flux function 

2 2 '=2 2 

where V^'f denotes the reconstructions of the conservative variables V from left 

and right sides of each cell interface (the reconstruction procedure employed in the 
present study will be described below). Consequently, the semi-discrete scheme takes 
the form 



dV, 

+ 



-^[^..k-^^-^] = (5-4) 



dt Ax' 

In order to discretise the advective flux F(V), we follow the method of [25, 26] and 
use the following FVCF scheme 

F(V) + F(W) , , F(W) - F(V) 
F(V,W) = ^ — - - U(V,W)-^ — ^— — 

The first part of the numerical flux is centered, the second part is the upwinding 
introduced through the Jacobian sign-matrix I[J(V, W) defined by 

U(V,W) = sign[A(i(V + W))], sign(A) = • diag(s+, s") • L, 

where = sign(A*). After some simple algebraic computations, one can find 

1/ S+ + S- {H/cs)is+ - s-)\ 



U 



2UcjH){s^-s-) 



the sign-matrix U being evaluated at the average state of left and right values. 

Finally the source term Sb{x,t) - (0, ^hh^tt), which is due to the moving bottom, is 
discretised by evaluating the bathymetry function and its derivatives at cell centers: 



^ ^b{x,t)dx [O, lh{xi,t) h^ttixi,t)) . 



A. 

Recall that the bathymetry is composed of the static part and of the landshde subject 
to a translational motion: 

h{x,t) = ho{x) - C(x,t) = ho{x) - Co{x -Xc(t)). 
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The derivative hxtt can be readily obtained from the formula 

5.1. High-order reconstruction. In order to obtain a higher-order scheme in space, 
we need to replace the piecewise constant data by a piecewise polynomial represen- 
tation. This goal is achieved by various so-called reconstruction procedures such as 
MUSCL TVD [36, 52, 53], UNO [34], ENO [33], WENO [57] and many others. In 
recent studies on unidirectional wave models [19] and on Boussinesq-type equations 
[18], the UN02 scheme showed a good performance with small dissipation in realistic 
propagation and run-up simulations. Consequently, we retain this scheme for the 
discretisation of the advective flux of the Peregrine system (5.1), (5.2). 

The main idea of the UN02 scheme is to construct a non-oscillatory piecewise- 
parabolic interpolant Q(x) to a piecewise smooth function V(x) (see [34] for more 
details). On each segment containing the face x^^i e [xi,Xi+i], the function Q(x) = 
is locally a quadratic polynomial and wherever v{x) is smooth we have 

Q(x) - V(x) = + 0(Ax3), ^(x±0) - ^ = + O(Ax^). 

ax ax 

Also, Q(x) should be non-oscillatory in the sense that the number of its local extrema 
does not exceed that of V(x). Since qj+i(a;j) = Vj and q^^i (xj+i) = Vj+i, it can be 
written in the form 

qi^ix) = V, + c),,i{V}x^— + iD,,i{V}x ^-^ , 

where Oj+i {V} = Vj+i - Vj and V is closely related to the second derivative of 
the interpolant since D^^i{V} = Ax^q" The polynomial c[-^i(x) is chosen to 

be the least oscillatory between two candidates interpolating V(a;) at (xj-i, Xj, Xj+i) 
and {xi,Xi+i,Xi+2)- This requirement leads to the following choice of S)j^i{V} = 

minmod(Di{V},Di+i{V}) with 

D.{V} = V,+i - 2V, + V,_i, D.+i{V} = V,+2 - 2V,+i + Vi, 

and where minmod(x,?/) is the usual minmod function defined as 

minmod(x,iy) = | [ sign(x) + sign(y) ] x min(|x|, 

To achieve the second order O^Ax"^) accuracy, it is sufficient to consider piecewise 
linear reconstructions in each cell. Let L{x) denote this approximately reconstructed 
function which can be written in this form 

L(x) = V, + S, xe[xii,Xi^]. 

/\X 2 2 

In order to L{x) be a non-oscillatory approximation, we use the parabolic interpola- 
tion Q(x) constructed below to estimate the slopes Sj within each cell 

Si = Ax X minmodf — ^ (xj - 0) , — ^ (xj + 0) ) . 

^ dx dx 
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In other words, the solution is reconstructed on the cells while the solution gradient 
is estimated on the dual mesh as it is often performed in more modern schemes [4, 5]. 
A brief summary of the UN02 reconstruction can be also found in [18, 19]. 

5.2. Treatment of the dispersive terms. In this section, we explain how we treat 
numerically the dispersive terms of the Peregrine system (5.1), (5.2) which are present 
only in the momentum conservation equation (5.2). We propose the following approx- 
imation for the second component of M(V) of M(V). 

M (\r\ - -h ~ '^hi{ut)i + hi_i{ut)i-i _ 1 (^f )»+i ~ 2(^t)» + (^f 

hj /t lT^.^ /- 2-o\,, hj (- 1- 



Note that this spatial discretisation is of the second order 0{Ax^) so as to be con- 
sistent with the UN02 advective flux discretisation presented above. If we denote by 
I the identity matrix, we can now rewrite the semi-discrete scheme in the form 

^.^[F«(V)-Fi^)(V)] = 0, 



where F± ' (V) are the two components of the advective numerical flux vector F at 



(2'\ 

the right (+) and left (-) faces correspondingly, and SJ^ denotes the discretisation of 



the second component of S;,. 

5.3. Time stepping. We assume that the linear system of equations is already in- 
verted and we have the following system of ODEs: 

Vi=AA(V,t), V(0)=Vo. 

In order to solve numerically the last system of equations, we apply the Bogacki- 
Shampine method proposed by Przemyslaw Bogacki and Lawrence F. Shampine in 
1989 [8]. It is a Runge-Kutta scheme of the third order with four stages. It has an 
embedded second order method which is used to estimate the local error and thus, 
to adapt the time step size. Moreover, the Bogacki-Shampine method enjoys the 
First Same As Last (FSAL) property so that it needs approximately three function 
evaluations per step. This method is also implemented in the ode23 function in 
Matlab [46]. The one step of the Bogacki-Shampine method is given by: 

= AA(V(") + |At„A;i,t„ + |At), 

k, = AA(V(")) + |At„A;2,t„ + |At), 

y(n+i) ^ v(")+At4|fci + iA;2 + |A;3), 

k^ = AA(V("+i),t„ + At„), 

V("+1) = v(")+At4^A;i + iA;2 + |A;3 + iA;4). 
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Convergence rate in norm 



Convergence rate in L2 norm 



—9 — FV scheme, slope = 1 .995 
2nd order 



- FV scheme, slope - 2.446 

- 2nd order 



Figure 6. This figure shows the convergence rate of the finite- volume 
scheme in the L°°-norm (left panel) and the L^-norm (right panel). The 
numerical integration of a solitary wave as shown in Figure 4 is compared 
to a translated profile. It appears that the second-order convergence is 
achieved. 

Here V(") V(t„), At is the time step and V2"'^^'' is a second order approximation to 
the solution V(t„+i), so the difference between V("+^) and V2"^^'' gives an estimation 
of the local error. The FSAL property consists in the fact that is equal to ki in 
the next time step, thus saving one function evaluation. 

If the new time step Atn+i is given by Atn+i - Pn^tn, then according to the H211b 
digital filter approach [47, 48], the proportionality factor p„ is given by 

Pn = [-) i-^) Pn% (5.5) 

where e„ is a local error estimation at time step and the constants [32 and a 
are defined by 

The parameter p gives the order of the scheme, and p = 3 in our case. 

Remark 1. The adaptive strategy (5.5) can he further improved if we smooth the 
factor pn before computing the next time step Atn+i'- 

At„+i = /3„At„, Pn = w(p„). 

The function oj{p) is called the time step limiter and should be smooth, monotonically 
increasing and should satisfy the following conditions: 

a;(0)<l, tj(+oo)>l, = = 1. 

One possible choice was suggested in [4:8]: 

uj(p) = 1 + /tarctanf^ ). 

In our computations the parameter k is set to 1. 
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5.4. Validation. The scheme described in this section is implemented in MATLAB, 
and runs on a workstation. To check whether the implementation is correct, we 
use the approximate solitary waves of (3.1), computed in the last section. These 
are used as initial data in the fully discrete scheme, and integrated forward in time. 
The computed solutions are then compared to the same solitary waves, but shifted 
forward in space by cto, where c is the wave speed, and to is the final time. This 
procedure is repeated a number of times with different spatial gridsizes. As a result, 
it is possible to find the spatial convergence rate of the scheme. As is visible in Figure 
6, the convergence achieved by the practical implementation of the discretisation is 
very close to the theoretical convergence rate. Since the temporal discretisation is 
adaptive, we do not present a convergence study in terms of the timestep. 



Let us consider a one-dimensional computational domain / = [a,b] = [0,220] com- 
posed of two regions: the generation region and a sloping beach on the right. More 
specifically, the static bathymetry function hQ(x) is given by a smoothed out profile 
generated from the expression 



In essence, this function represents a perturbation of the sloping bottom by two un- 
derwater bumps. We made this nontrivial choice in order to illustrate the advantages 
of our landslide model, which was designed to handle general non-fiat bathymetries. 
The parameters can be chosen in order to fit a given bathymetry, but the particular 
values used here are Ai = 4.75, A2 = 8.85, ki = 0.06, ^2 = 0.13, xi = 45, X2 = 80, and 
m = 120. The bottom profile for these parameters is depicted in Figure 7. Of course, 
in general, if the bottom topography is known, then a numerical bathymetry map 
could also be used. 

We now present some results of the solution of the surface wave problem using the 
model in Section 3, integrated numerically with the method of Section 5. A landslide 
is introduced on the left side of the bathymetry, and using the method of Section 2, 
its path along the bottom is determined by following the barycenter. Simultaneously, 
the system (3.2) is solved with the time- dependent bottom topography given from 
the solution of the landslide problem. The problem is integrated up to a final time 
T. Figure 8 shows wave records at six virtual wave gauges for both the dispersive 
system (3.2) and the shallow-water system. It appears from this figure that the 
shallow-water system underpredicts the development of free-surface oscillations. In 
particular, the wave gauges located at a; = 40 and x = 60 show similar waveheights for 
both the shallow-water, and the dispersive system, but a qualitative divergence, as 
small oscillations are already developing which are not captured by the shallow-water 
system. Once the waves have propagated to the wave gauges located at x = 80, the 



6. Numerical results and discussion 



[ c?o + tan 5 • (m - 
where the function p(x) is defined as 
p(x) = Aisech {ki{x 




I) + p{x), a < X <m 

a) - tan d ■ {x - m) , x > m, 



Xi)) + A2sech {k2(x - X2)) 
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Symbol 


Parameter 


Units 


Values 


a 


gravitational acceleration 


m/s 


9.81 


(in 


water depth at x = a 


m 


1.0-2.0 


tan((5) 


bottom slope 




0.1 


A 


landslide amplitude 


m 


0.55 


I 


landslide length 


m 


52.4 




added mass coefficient 




1.0 


Cd 


water drag coefficient 




1.0 


c/ 


friction coefficient 




tan(3°) 


7 


density ratio water /landslide 




1.8 




friction coefficient with bottom 




7.63e-4 


Cy 


viscous friction coefficient 




1.27e-3 



Table 1. Values of various parameters used in numerical computations. 




Figure 7. This figure shows the physical setup of the problem. The river 
bed is indicated in dark grey. The computational fluid domain is shaded light 
grey, and the landslide is visible in black. Note the difference in horizontal 
and vertical scales in the left panel. The right panel shows a closeup of the 
left beach and part of the landslide in a one-to-one aspect ratio. 



dispersive oscillations have amplified, so that the waveheight is larger by a factor of 
2 to 3 than the waveheight predicted by the shallow-water system. Going further to 
the wave gauges located at x = 100 and x = 120, the now rising bottom starts to have 
a damping effect on the waves. The maximum and minimum free-surface elevation 
is shown in Figure 9. 

Figure 10 shows the development of the kinetic energy of the landslide mass and 
simultaneously the total (kinetic plus potential) energy contained in the body of the 
fluid and the surface waves. Energy development is an important question in the 
study of tsunamis, and there have been studies exclusively devoted to this question 
[50]. Energy issues connected to water wave models of Boussinesq type have also 
been studied before [1, 2, 16]. While these models contained a source of energy, in 
the case at hand, the work done by friction as the landslide slides down the bottom 
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Figure 8. This figure shows time series of the surface elevation at wave 
gauges located at x - 40, x - 60, x - 80, x - 100 and x - 120. The solid (blue) 
curve depicts the wave elevation computed with the dispersive system (3.2), 
and the dashed curve represents results obtained from the shallow-water 
system. All variables are non-dimensional. 



acts as a drain of energy, and after the landslide has come to rest, all energy has been 
transferred to the fluid. However, not all energy can be considered as residing in the 
wave motion, because a significant amount of energy is needed to lift the water from 
the final position of the landslide to the initial position of the landslide. This results 
in a large increase in potential energy of the fluid, and only a fraction of the potential 
energy of the landslide is transferred to the wave motion. This fact has also been 
explained in previous works [32]. 
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Maximum absolute value el the tiorizontal velocity 
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Figure 9. This figure shows the maximum and minimum of the surface 
excursion, and the horizontal velocity as a function of (non-dimensional) 
time. 



In order to compute the wave energy in the fluid, we use the integral 

= £{^V^ + liho + VW] dx, (6.1) 

which arises from the shallow-water theory. The kinetic energy of the landslide is 
given by 

Esi = \mv\ (6.2) 

with the generalised mass m given by (2.3), and f = ^ as deflned in Section 2. Figure 
10 shows the development of the wave energy and kinetic energy of the landslide. The 
upper panel shows the energy according to the shallow-water and dispersive model. 
The lower panel shows the kinetic energy of the landslide. We have also computed 
the Froude number Fr = , " ^ during the evolution. Here v is the x-component of 

the velocity of the barycenter of the landslide, Xc is the position of the barycenter, 
and h{xc) is the corresponding local water depth. This number was always found to 
be much less than 1 in all numerical experiments. The maximum value was generally 
about 0.5. To compute the wave run-up and draw-down, we use exact representations 
given by Choi et. al. [11] (a similar formula was also derived in [15]). On the right 
beach, the undisturbed water depth at the edge of the computational domain is 
h - 3, and the distance from the computational domain to the shore line is L = 30. 
Using the shallow-water wave speed, the travel time of a wave from the edge of the 
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Water wave energy evolution 
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Figure 10. This figure shows the development of the wave energy, and the 
kinetic energy of the landshde as a function of (non-dimensional) time. Note 
that the kinetic energy of the landslide starts from (all energy is potential) 
and also ends at (all energy has been dissipated or transferred to the fluid). 

Run-up on the left beach 



Peregrine 




50 100 150 200 250 



Run-up on the right beach 

1 ^ \ \ 




Figure 1 1 . This figure shows the run-up on the left and right beach using 
(6.4), computed with the dispersive system (solid curve) and the nonlinear 
shallow-water system (dashed curve) as a function of (non-dimensional) time. 
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computational domain to the shore is computed as 



Vgh V get 

Then the formula for the wave run-up R at the shore reads 




(6.3) 



Jo 



t-T 



t-T drj 



(t - r)2 - T2 dT 



(x, t) dT 



(6.4) 



with X = 220. At the left beach, the undisturbed water depth is h = 1.642, and the 
distance to the beach is L = 11.2814. A similar formula can be then be computed for 
X = 0. 

Figure 11 shows the run-up on the left and right beaches both in the Boussinesq 
scaling and in the shallow-water theory. While the agreement is fair on the left beach, 
it appears immediately that the Boussinesq theory predicts a wave run-up on the right 
beach which is much larger (roughly by a factor of two) than the wave run-up ac- 
cording to the shallow-water theory. A possible explanation for this divergence is 
the nature of the numerical solver when applied to the shallow-water system. In this 
case, there is continuous numerical dissipation through the handling of hyperbolic 
wave breaking. Since the waves do not break in the Boussinesq scaling, the dissi- 
pation is not present, or at least much smaller. The difference can also be read off 
from the comparison of the wave energy in the Boussinesq and shallow-water system 
provided in Figure 10. It can be seen there that the wave energy in the shallow-water 
model starts to diverge from the Boussinesq model at non-dimensional time t = 50. 
The difference between the two increases continuously, until at the final time, the 
Boussinesq energy is about 50% larger than the shallow-water energy. Note that 
significant run-up in Figure 11 does not happen until non-dimensional time t - 75, 
at which time the energy in the Boussinesq system is already much larger than in 
the shallow- water system. In Figure 12, we have plotted the maximum wave am- 
plitude, the minimum wave amplitude, and the maximum wave run-up on the left 
and right beaches. In comparison to previous studies, such as [30], where an open 
domain was used, it appears that in our case, the maximal amplitude, as well as the 
run-up have a minimum at do between 1 and 1.5. In [30], it was found that maximum 
wave amplitude and run-up (on the left beach) were strictly decreasing functions of 
do- The phenomenon of rising amplitude and run-up may be accredited to resonant 
effects which are absent on an open domain (such as an ocean beach), but cannot be 
neglected for tsunamis generated by landslides in rivers and lakes. 



The influence of an underwater landslide on surface waves in a closed basin have 
been studied. The key features of the study have been that the motion of the underwa- 
ter landslide have been determined by integrating a second-order ordinary differential 
equation derived from first principles of Newtonian mechanics, and that the wave mo- 
tion has been studied in the Boussinesq scaling which allows for both nonlinear and 
dispersive effects. The dynamics of the motion of the bottom have been developed 



7. Conclusion 
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Figure 12. This figure shows the maximal and minimal wave amplitude, 
and the maximum run-up on the left and right beaches as a function of the 
initial depth of the center of the landslide do . 

following recent work in [35]. The Boussinesq model which has been utilised here 
allows for a dynamic bathymetry, and was derived in [56]. The numerical method 
used in this paper is an extension of the method put forward in [4, 5]. The results 
presented in Section 6 clearly show that dispersion may have a strong effect on the 
run-up and draw-down at the beaches. Of course, this difference could be more or 
less pronounced depending on the particular case under study. For example, the di- 
vergence between the shallow-water theory and the dispersive model is stronger at 
the right beach than at the left beach. The results also show that a finite domain 
exhibits different behavior than a half-open domain (such as used in [30]) with respect 
to the dependence of the wave run-up on the initial depth of the landslide. While 
the run-up is a strictly decreasing function of the initial depth in an open domain, a 
closed domain appears to exhibit resonant effects, which make the dependence more 
complex. 
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